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Abstract 

The energy levels of neutral atoms supported by Yukawa potential, 
V(r) = — Zexp(— ar)/r, are studied, using both dimensional and dimen- 
sionless quantities, via a new analytical methodical proposal ( devised 
to solve for nonexactly solvable Schrodinger equation). Using dimension- 
less quantities, by scaling the radial Hamiltonian through y = Zr and 
a = a/Z, we report that the scaled screening parameter a is restricted to 
have values ranging from zero to less than 0.4. On the other hand, working 
with the scaled Hamiltonian enhances the accuracy and extremely speeds 
up the convergence of the energy eigenvalues. The energy levels of several 
new eligible scaled screening parameter a values are also reported. 



1 



1 Introduction 



The Yukawa potential plays special roles in different branches of physics. It is 
known as the Dcbyc-Hiickel potential in plasma physics, as the Thomas-Fermi or 
screened Coulomb potential in solid - state and atomic physics, and the Yukawa 
potential in nuclear physics. It is also used in the models describing metal - in- 
sulator transition [1,2]. As such, several approaches have been invested to find 
the energy levels of this potential. Amongst exist the simple variational proce- 
dure [3] , the use of atomic orbitals with a set of fitting parameters [4] , Rayleigh 
- Schrodinger perturbation theory [5,6], method of potential envelopes [7], iter- 
ative procedure [8], different numerical calculations [9,10], and the shifted 1/N 
expansion method [11] based on Papp's [12] or Imbo et al.'s [13,14] proposals. 
On the other hand, the Yukawa form 

V(r) = —Zexp(—ar)/r (1) 



is known to simulate the screening potential of neutral atoms, where the screen- 
ing parameter a is chosen to be a = a. Z x ^ corresponding to the Z - depen- 
dence of the reciprocal of the Thomas-Fermi radius of the atom. However, 
invoking Fermi-Amaldi correction [15] in the context of Ecker-Weizel approxi- 
mation (EWA) method [16], Dutt and Varshni [17] have amended the screening 
parameter and suggested that a — a Z x ^(\ — 1/Z) 2 ^ 3 with a a = 0.98. The 
correctness of this amendment has been justified by Lai and Madan [18], who 
have used the hypervirial - Pade method (HPM) to obtain very accurate energy 
eigenvalues. Therefore, the energy levels of the neutral atoms, supported by 
potential (1), have been reexamined by applying the shifted large - N expansion 
(SLNT) ( where N is the number of spatial dimensions) [14]. However, some 
drawbacks of the HPM and SLNT are in order. The HPM involves elaborate 
computational time and effort for each numerical prediction, and its application 
becomes quite restricted due to the absence of compact analytical expressions 
for the eigenenergies and eigenfunctions. And, because of the complexity in 
handling large - order corrections of the standard Rayleigh - Schrodinger per- 
turbation theory, only low - order cigenencrgy calculations have been reported 
for SLNT and large order - calculations have been neglected [11-14]. Yet, nci- 
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ther SLNT nor HPM is utilitarian in terms of calculating the eigenvalues and 
eigenfunctions in one batch. Moreover, Scherrer et al. [19] have concluded, via 
matrix continued fractions method, that SLNT leads to dubious accuracies. A 
conclusion that we have very recently confirmed [20] . 

In this work we introduce a new methodical proposal for solving Schrodingcr 
equation with an arbitrary spherically symmetric potential. We shall use l/l as 
an expansion parameter, where / = I — ft, I is a quantum number and (3 a suitable 
shift introduced, mainly, to a avoid the trivial case / = 0. In fact, the scope of 
its applicability is not limited to the spherically symmetric problems, it extends 
also to systems at lower dimensions. Here is the conversion strategy. One would 
rewrite the centrifugal term 1(1 + l)/2r 2 , in Eq.(2) below, as A(A + l)/(2q 2 ) and 
use 1/A as an expansion parameter. In this case, A = I in three dimensions, 
where I is the angular momentum quantum number with q = r > 0, A = 
|m| — 1/2 in two dimensions, where m is the magnetic quantum number with 
q = (x 2 +y 2 ) 1 ^ 2 > 0, and A = —1 and/or in one dimension with — oo < q < oo. 
Thus the method could be called a pseudoperturbative shifted - I expansion 
technique (PSLET), to be described in the following section. 

2 The Method 

To construct the method we start with the time - independent one - dimensional 
form of Schrodingcr equation, in % = m = 1 units, 



2dr 2 + ( 2r 2 ^ + 



*n P ,i(r)=^ P ,i*„ P ,i(r). (2) 



and to avoid the trivial case 1=0, the quantum number I is shifted through the 
relation 1 = 1 — fi. Eq.(2) thus becomes 



Id 2 - 1 



(3) 



- (r) = p + ( 2/?+ i)r + ^ + i) + P m (4) 
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Herein, it should be noted that Q is a constant that scales the potential V(r) 
at large - I limit and is set, for any specific choice of I and n r , equal to P at the 
end of the calculations [11-14,20-22]. And, [3 is to be determined in the sequel. 

Our systematic procedure begins with shifting the origin of the coordinate 
through 

.T = F/ 2 (r-r )/r , (5) 



where r$ is currently an arbitrary point to perform Taylor expansions about, 
with its particular value to be determined. Expansions about this point, x = 
(i.e. r = r ), yield 



1 

12 ~~ 



>+l) 



n=0 



x n l 



-n/2 



(6) 



V(x(r)) = 



n=0 



n V(r ) \ (r g)" -_ w/2 
rfr^ y n! 



(7) 



Obviously, the expansions in (6) and (7) center the problem at an arbitrary 
point r and the derivatives, in effect, contain information not only at r but 
also at any point on the axis, in accordance with Taylor's theorem. Also, it 
should be mentioned here that the scaled coordinate, equation (5), has no effect 
on the energy eigenvalues, which are coordinate - independent. It just facili- 
tates the calculations of both the energy eigenvalues and eigenfunctions. It is 
also convenient to expand E as 



(8) 



Equation (3) thus becomes 



■5F + f^» 



*n P ,lO»0 = -jE„ r ,l^>„ r ,l{x), 



(9) 
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with 



■V(x(r)) - r 2 J 



1 . V(ro) 
H Q 



+ p/ 2 



-x + 



V (r )r%x 

Q 



+ 



3 2 V"(r )rW 
2 + 



2Q 



+ ( 2/ 9 +!)£(-!)' 



(n+1) 



n=l 



+ ri 



°^jL ( } 2rg U dr™ ) nlQ 



J-(n-2)/2 



+ w + i)f](-i)«i^ a; "r-(™+ 2 )/ 2 + ^±^, (io) 



n=0 



where the prime of V(ro) denotes derivative with respect to ro- Equation (9) is 
exactly of the type of Schrodinger equation for one - dimensional anharmonic 
oscillator 



\^ + \ w ^+e + P{x) 



X„ r (x) = X nr X n Jx), 



(11) 



where P(x) is a perturbation - like term and eo is a constant. A simple com- 
parison between Eqs.(9), (10) and (11) implies 



£o = I 



1 , rgy(rp) 

2 Q 



2/3+1 0(0 + 1) 
2 + 2F ; 



(12) 



and 



1 

+ I 



"1 . rgy(rp) 
2 + Q 

/?(/? + !) 



+ 



2/9+1 , , 1. 
— — + K + 



71 = 2 



(13) 



n=-2 



(14) 
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Equations (13) and (14) yield 



( _ 2) = i 

2r 2 Q 



(15) 



2/3+1 , , , \ ' 



(16) 



£ (0) = 1 



/?(/? + !) 



(17) 



(18) 



f — 2) 

Here ro is chosen to minimize £^ / , 



i. e. 



dr 



= and 



■2 P (-2) 



drl 



> 0. 



(19) 



Hereby, V(r) is assumed to be well behaved so that 2?( 2 ) has a minimum ro and 
there are well - defined bound - states. Equation (19) in turn gives, with I = ^Q, 



1-13= yJr*V'(r Q ). 



(20) 



Consequently, the second term in Eq.(10) vanishes and the first term adds a 
constant to the energy eigenvalues. It should be noted that energy term PE^ 2 ^ 
has its counterpart in classical mechanics. It corresponds roughly to the energy 
of a classical particle with angular momentum L z =l executing circular motion 
of radius r in the potential V(ro). This term thus identifies the leading - order 
approximation, to all eigenvalues, as a classical approximation and the higher 
- order corrections as quantum fluctuations around the minimum ro, organized 
in inverse powers of I. 
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The next leading correction to the energy series, lE^W consists of a con- 
stant term and the exact eigenvalues of the unperturbed harmonic oscillator po- 
tential w 2 x 2 /2. The shifting parameter (3 is determined by choosing F^~^=0. 
This choice is physically motivated. It requires not only the agreements between 
PSLET eigenvalues and the exact known ones for the harmonic oscillator and 
Coulomb potentials but also between the eigenfunctions as well. Hence 



P = - 



1 / ^ 

- + K + -)w 



(21) 



where 



w = \l 3 + 



r V"(r ) 



V'(r ) 



(22) 



Then equation (10) reduces to 



] -V(x(r)) = rll 



1 , V(ro) 
H Q 



+ £y")(z)r 



■n/2 



n=0 



(23) 



where 



(0) I \ 1 2 2,2/3+1 
v ( >(x) = -w £ x z + - , 



(24) 



vW(x) =-(2(3+l)x- 2x 3 + ^ M ~ 3 



6Q 



(25) 



and for n > 2 



vW(x) = (-ir(2/3 + l)^±^x" + (-l)«^±ll(n-l) a ;("- 2 ) 



+ 



, lV > + 3) r<" +4) d^V(r ) 
[ ' 2 + Q(n + 2)! dr-+ 2 



r.n+2 



(26) 
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Equation (9) thus becomes 



n=0 



*»r,lO»0 = 

\ ' n— 9 



*»r,lO»0- 



(27) 



When setting the nodclcss, n r = 0, wave functions as 

*o,i(x(r)) = exp(U ,i(x)), (28) 
equation (27) is readily transformed into the following Riccati equation: 

- ^[u"(x) + u\x)u\x)\ + J2 v(n) ( x )i~ n/2 

n=0 



Hereinafter, we shall use U(x) instead of Uq^(x) for simplicity, and the prime 
of U(x) denotes derivative with respect to x. It is evident that this equation 
admits solution of the form 

oo oo 

U'(x) = J2 U (n) (x)r n / 2 + G {n \x)l' (n+1)/2 , (30) 

n=0 n=0 

where 

n+l 

U^(x) = Y. D ^ x2m ^ > D °,n = 0, (31) 
m=0 

n+l 

G^(x)=J2 C m, n x 2m . ( 32 ) 

m=0 



1 {(3 (13 + 1) (o) 
— o !" A 



+ H A o l~ n - ( 29 ) 



n=2 
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Substituting equations (30) - (32) into equation (29) implies 



_ 1 ^ \^j(n)' J-n/2 + Q(n)' J-(n+l), 

^ n=0 

oo oo 

_ -^2^2 \u {n) u {p) i- {n+p)/2 + G ( - n) G {p) r {n+p+2)/2 + 2U { - n) G {p) r {n+p+1)/ '- 1 



n=0p=0 

+ ±v^n-^^l^^^ + x^y±xr i} r n , (33) 

IT, — IT, — 2 

where primes of U^ n \x) and G^ n \x) denote derivatives with respect to x. 
Equating the coefficients of the same powers of / and x, respectively, ( of course 
the other way around would work equally well) one obtains 

_i^(o)'_ ^(0)^(0) +v (o )=0) (34) 



U<-°) (x) = Di,o ; Di,o - -w, (35) 
and integration over dx yields 

XjW{x) = -wx. (36) 

Similarly, 

_ ![[/(!)' + G (o)'j _ u(°)uW _ [/(o) G (o) + = 0j (37) 

U w (x) = Q, (38) 

G (0) (^) = Co,o + C 1 , 2; 2 , (39) 
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(40) 




10 



(49) 



• • • and so on. Thus, one can calculate the energy eigenvalue and the eigen- 
functions from the knowledge of C mj „ and D m n in a hierarchical manner. Nev- 
ertheless, the procedure just described is suitable for systematic calculations 
using software packages (such as MATHEMATICA, MAPLE, or REDUCE) to 
determine the energy eigenvalue and eigenfunction corrections up to any order 
of the pscudopcrturbation series. 

Although the energy series, Eq.(8), could appear divergent, or, at best, 
asymptotic for small I, one can still calculate the eigenenergies to a very good 
accuracy by forming the sophisticated Pade' approximation to the energy series. 
Our strategy is therefore clear. 

3 Yukawa Potential 

Let us begin with the Yukawa potential (1), where the screening parameter a 
is given by the amended Z-dependent relation 

a = Q.98Z 1/3 (l-l/Z) 2/3 . (50) 



Thus, Eq.(22) with n r = reads 



w= / -^o+^o + l 

V ar + l y ' 



and Eq.(20) in turn implies 



I + ^(1 + w) = ^Zr e- ar °(ar + 1). (52) 



Evidently, one should appeal to numerical techniques to solve for r since a 
closed form solution for Eq.(52) is hard to find, if it is feasible at all. Once r 
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is determined the coefficients C m ,n and D m . n are obtained in a sequel manner. 
Consequently, the eigenvalues and eigengfunctions are calculated in the same 
batch for each value of Z and a. 

In tables 1 and 2, we list our K- and L-shell binding energies, measured in 
units of 2Ry=27.196 eV, for some values of Z and compare with those of the 
hypervirial - Pade method (HPM) [18], and of the shifted large - N expansion 
technique (SLNT) [14]. In table 1, the results of Ecker - Weizel approximation 
(EWA) [16] are also displayed. The first ten terms of our energy series are 
calculated and the stability of the sequence of Pade approximants is considered. 

Evidently, our results compare better with those of HPM than the results of 
SLNT. However, considering only the first four terms of our series, our results 
are found in exact agreement with those of SLNT. The accuracy of our results 
is enhanced by the Pade approximants which show stability all the way from 
E[2, 2] to E[4, 5]. Such stability indicates that our results are exact, provided 
that some uncertainty lies in the last j digits in parentheses. 

On the other hand, one would pass to dimensionless quantities by scaling 
the Hamiltonian and reducing it to the form 

where y G (0, oo) and a 1 is a positive screening parameter. Then, to return back 
to dimensional quantities one has to make the following substitutions 

V = Zr , a ' = !=0.98Z- 2 / 3 (l-l) 2 / 3 , (54) 

where Z — 1, 2, 3, • • •. Therefore, a 1 can not be a random screening parameter. 
It is restricted by Eq.(54), which in effect implies that the lowest value of a 1 is 
zero when Z = 1 and immediately a maximum value 0.38891326 at Z—2. For 
Z > 2, a' decreases smoothly and converges to zero as Z — > oo. To figure this 
out, one should simply sketch the relation of a' vs Z, given in Eq.(54). Hence, 
we proceed with our comparison, using only several eligible values for a', and 
list our results (tables 3 and 4) along with those of Varshni [11], via SLNT, 
Fernandez et al. [23], via a rational function approximation (RFA), and the 
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exact ones [24,25]. Our results, again, are in better agreement with the RFA 
and the exact ones than those of Varshni [11], via Papp's [12] or Imbo's [13] 
SLNT proposals. 

Although we have used the atomic units through out, our energies in tables 
1 and 2 are multiplied by 2Ry for the sake of comparison with the results from 
HPM [18], EWA [16], and SLNT [14]. 

4 Concluding remarks 

In view of the results reported in tables 1-4, some remarks merit to be mentioned. 

It is evident that working with the scaled Hamiltonian , Eq.(53), facilitates 
the calculations and enhances the accuracy of PSLET. As can be seen from our 
results, the convergence for the energy eigenvalues is extremely rapid when the 
scaled Hamiltonian is used. It should be obvious for the reader to see that all 
Z - values in tables 1-2 fall in the range of a' in tables 3 and 4. Moreover, as 
a' gets smaller the Yukawa potential becomes more Coulombic in nature, and 
when a'=0 at Z=l the Yukawa collapses into the exact Coulombic potential, for 
which PSLET gives exact spectrum from the leading term pE^-~ 2 ^ and higher - 
order corrections vanish identically. 

The answer to the question as to why all methods reported by Varshni [11] 
give poor accuracies for a' > 0.7, lies perhaps in the restriction provided by 
Eq.(54) that a' < 0.4 (table 1 of ref.[ll] bears this out). Moreover, the non- 
negligible deviations of SLNT results [19,20,26-30] from the exact values should 
be attributed, mainly, to the limited capability of SLNT to calculate the energy 
corrections beyond the fourth - order. 

From the knowledge of C m . n and D mn one can calculate, in the same batch, 
the wave functions to study electronic transitions and multiphoton emission 
occurring in atomic systems in the presence of intense laser fields, for example. 
Such studies already lie beyond the scope of our present methodical proposal. 

To sum up, we have demonstrated that it is an easy task to implement 
PSLET without having to worry about the ranges of couplings and forms of 
perturbations in the potential involved. In contrast to the textbook Raylcigh - 
Schrodinger perturbation theory, an easy feasibility of computation of the eigen- 
values and eigenfunctions, in one batch, has been demonstrated, and satisfactory 
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accuracies have been obtained. Moreover, a nice numerical trend of convergence 
has been achieved. Nevertheless, another suitable criterion for choosing the 
value of the shift /?, reported in ref. [30] , is also feasible. This reference should 
be consulted for more details. 
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Table 1: K-shell energies, with n r — and I = 0, in KeV for V(r) = —Ze~ ar /r. 
Where Epslet represents PSLET results, Ehpm denotes the Hypervirial 
method [18], Eewa from [16], and Eslnt from [14]. £[4,5] is the [4,5] Pade 
approximant obtained by replacing the last j digits of £[4,4] with the j digits 
in parentheses. 



z 


—Eewa 


—Ehpm 


—Eslnt 


— E PS LET 


-£[4,4] & (-£[4,5]) 


3 


0.05334 


0.05415 


0.05414 


0.054135 


0.05414708 (19) 


6 


0.274 


0.27623 


0.27623 


0.276226 


0.27622917 (22) 


9 


0.701 


0.70461 


0.70461 


0.704612 


0.70461288 (9) 


14 


1.897 


1.90320 


1.90320 


1.903198 


1.903198386 (6) 


19 


3.716 


3.72545 


3.72545 


3.725447 


3.725447269 (7) 


24 


6.171 


6.18277 


6.18277 


6.182765 


6.182765210 (07) 


29 


9.268 


9.28212 


9.28213 


9.282124 


9.282124421 (18) 


34 


13.012 


13.02830 


13.02830 


13.028304 


13.028303677 (4) 


39 


17.407 


17.42482 


17.42482 


17.424818 


17.424817591 (88) 


44 


22.454 


22.47438 


22.47438 


22.474378 


22.474378305 (3) 


49 


28.157 


28.17915 


28.17915 


28.179153 


28.179153161 (0) 


54 


34.517 


34.54092 


34.54092 


34.540921 


34.540920826 (5) 


59 


41.535 


41.56612 


41.56117 


41.56117189 


41.561171872 (1) 


64 


49.213 


49.24118 


49.24118 


49.2411768 


49.241176745 (4) 


69 


57.553 


57.58203 


57.58203 


57.5820335 


57.58203349 (8) 


74 


66.554 


66.58470 


66.58470 


66.5847023 


66.584702299 (9) 


79 


76.217 


76.25003 


76.25003 


76.2500313 


76.250031268 (8) 


84 


86.544 


86.57878 


86.57878 


86.5787759 


86.578775902 (1) 
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Table 2: L-shell energies, with n r — and I = 1, in KcV for V(r) = — Ze~ ar jr. 
Where Epslet represents PSLET results, Ehpm denotes the Hypervirial 
method [18], and Eslnt from [14]. E[A, 5] is the [4,5] Pade approximant ob- 
tained by replacing the last j digits of E[A, 4] with the j digits in parentheses. 



z 


—Ehpm 


—Eslnt 


— EpsLET 


-£[4,4] & (-£[4,5]) 


9 


0.00423 


0.00437 


0.0056 


0.0041 (1) 


14 


0.08953 


0.08945 


0.08948 


0.089523 (31) 


19 


0.28398 


0.28394 


0.283962 


0.283977 (8) 


24 


0.59991 


0.59989 


0.599905 


0.599912 (2) 


29 


1.04510 


1.04508 


0.045092 


1.0450957 (8) 


34 


1.62485 


1.62484 


1.624851 


1.6248535 (6) 


39 


2.34309 


2.34308 


2.343087 


2.3430888 (9) 


44 


3.20280 


3.20280 


3.203800 


3.2028012 (3) 


49 


4.20638 


4.20637 


4.2063780 


4.20637855 (9) 


54 


5.35577 


5.35577 


5.3557723 


5.35577266 (9) 


59 


6.65261 


6.65261 


6.6526133 


6.65261354 (6) 


64 


8.09829 


8.09828 


8.0982857 


8.09828593 (5) 


69 


9.69398 


9.69398 


9.6939829 


9.69398303 (4) 


74 


11.44075 


11.44074 


11.4407452 


11.44074536 (7) 


79 


13.33949 


13.33949 


13.3394894 


13.33948956 (7) 


84 


15.39103 


15.39103 


15.3910302 


15.39103029 (30) 



18 



Tabic 3: K-shell energies, with n r = and 1 = 0, inh = m = e = l 
units for V(r) = —e~ a r jr. Where Epslet represents PSLET results, Ep app 
denotes SLNT, via Papps proposal [11], and Ej m b SLNT, via Imbo's proposal 
[11]. £[4, 5] is the [4,5] Pade approximant obtained by replacing the last j digits 
of £[4,4] with the j digits in parentheses. 



a' 




-^Imbo 


— EpsLET 


-£[4,4] & (-£[4,5]) 


Exact [23-25] 


0.1 


0.407058 


0.407058 


0.40705803 


0.407058031 (1) 


0.407058031 


0.2 


0.326808 


0.326808 


0.326808231 


0.326808515 (9) 


0.326808515 


0.3 


0.257634 


0.257634 


0.257628168 


0.25763869 (8) 


0.25763869 


0.4 


0.198346 


0.198345 


0.198260722 


0.198377 (8) 


0.198377 



Table 4: K-shell energies, with n r = and Z = 0, in?i = m = e= l units for 
V(r) = -e- a ' r /r. Where E PSLET represents PSLET results, £[4, 5] is the [4,5] 
Pade approximant obtained by replacing the last j digits of £[4,4] with the j 
digits in parentheses. 



a 1 


— E PSLET 


-£[4,4] & (-£[4,5]) 


0.01 


0.4900745067 


0.490074506746694 (4) 


0.02 


0.48029610598 


0.48029610598378 (8) 


0.03 


0.4706620270266 


0.4706620270246 (4) 
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